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I. INTRODUCTION 

In this paper we study a variation of the one- 
dimensional Hubbard model (HM) in which the difference 
between the hopping amplitudes, say tf > t±, is respon- 
sible for an explicit breaking of the rotational symmetry. 
It is described by the Hamiltonian 



h.c. 



(1) 



where Cj a denotes the annihilation operator of a fermion 
with a =t, 4 a t site j and rij a = Cj a Cj a are the associated 
number operators. 

This asymmetric Hubbard model (AHM) has been 
studied in the past [f| to describe the essential features 
of the metal-insulator transition in rare-earth materials 
and transition-metal oxides; in this case a represents two 
types of spinless fermions (the real spin being considered 
not essential for the transition to be modelled). The 
"light" particles are described by band (Bloch) states, 
while the "heavy" ones tend to be localized on lattice 
(Wannier) sites. More recently, this model has gained a 
renewed interest in experiments with optical lattices, in 
which both the effective strengths of the kinetic and of 
the potential parts can be varied in a rather controlled 
way, including the possibility of reaching the attractive 
regime U < 0. The possibility to use cold atoms (13, [HJ 
to engineer condensed matter systems with a high tun- 
ability offers an experimental way to test theoretical re- 
sults with great accuracy. Two-species models with dif- 
ferent hopping coefficients can also be realized by trap- 
ping atomic clouds with two internal states of different 
angular momentum, thereby introducing a spin depen- 
dent optical lattice, which enables to modify the ratio 



a = ti/tf by controlling the depth of the optical lat- 
tice [![. Yet another possibility is to trap two different 
species of fermionic atoms, so that the "anisotropy" a is 
given naturally by the ratio of masses. In the context 
of cold atoms in optical lattices the subscripts a =f, 4- 
are not related to the electron spin but label the two 
different species of fermions, either different atoms with 
half-integer spin or different excited states of one atomic 
specie with fine structure splitting. 

Many recent papers on the subject are devoted to the 
onset of the so-called Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) phase, which may occur with unbalanced species 
[3, In this paper we will consider instead the case 
of balanced species: = N±. The parameters that 
influence the phase diagram can be cast in the form 
of an anisotropy coefficient z = (t-j- — f|,)/(i-f + t±) = 
(1 — a) /(I + a) and a dimensionless onsite potential 
u = U/t, where t = (t^ + tj,)/2 is an overall energy scale. 
In addition one can consider the effect of the total fill- 
ing n = N/L, N being the total number of fermions and 
L the chain length. At variance with the typical situa- 
tion in condensed matter physics, where the bulk filling 
and magnetization m z = {n^ — n^)/2 are controlled in 
a grand-canonical framework by the chemical potential 
and an external magnetic field respectively, in the con- 
text of cold atoms it is conceivable to fix independently 
the number of particles in each of the two species, de- 
spite the fact that only a given choice of the densities 
n-j- and might correspond to the absolute minimum of 
the grand potential. Henceforth we will assume that they 
are equally populated, that is = n\. = n/2 and we will 
limit ourselves to n < 1. Note that if E(n, m z ) is the en- 
ergy in a given sector of n and m z (these two quantities 
are good quantum numbers also in the asymmetric case), 
then a particle-hole transformation rij a — > (1 — nj a ) leads 
to E{2 - n, -m z ) = E(n, m z ) - UL(n - 1), so that for 
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m z = it is sufficient to analyze the phase diagram for 
positive and negative U and n < 1, in order to infer fea- 
tures also at n > 1. The phase diagram at half-filling, 
n = 1, has been studied in ref. [2l| , soon after the devel- 
opment of White's density-matrix renormalization group 
(DMRG) method @ . The limiting case of the symmetric 
HM (z = 0) can be solved exactly via the Bethe ansatz 
approach |9j . The opposite extremal case z = 1 is usually 
called the Falicov-Kimball (FK) model. In any dimen- 
sionality it has been proved (Iol | that for large positive U 
and the system has a ground state characterized by 
a spatially non-homogeneous density profile where the 
heavy particles are compressed in a region with length 
Li = N±. In this paper we refer to this phase as totally 
segregated state (TSS). A similar phase, where the heavy 
particles acquire a finite kinetic energy (Zq > N±), ap- 
pears for z < 1 0] and can be interpreted as the result of 
an effective attractive interaction between light fermions, 
mediated by the heavy ones (see [lij in ID and [l6j in 
2D). The latter situation is denoted in the following as 
phase separation (PS). A cartoon representation of the 
different phases is presented in Fig. [TJ At smaller asym- 
metries, the system is instead in a more conventional spa- 
tially homogeneous phase (HP). The ground state phase 
diagram for the model in the m z = sector has been 
discussed in [TJ by means of the bosonization approach. 
In this context [Tl| the HP-PS transition line at u > 
has been interpreted as the curve in (z,u)-plane where 
the velocity of one of the two decoupled bosonic modes 
vanishes. 

Let us summarize the content of this paper. In Sec. 
HIl we will study numerically the attractive regime (U < 
0), by using a DMRG program. In particular we will 
examine which kind of correlations (charge or pairing) is 
dominant in this region of the phase diagram. Then we 
will move to consider the repulsive (U > 0) regime. In 
Sec. IIHI we will analytically discuss the weak coupling 
regime (\U\ <C t a ) by means of a variational method that 
compares the energy of the PS state with that of the HP 
state, the latter being calculated within a second order 
perturbation analysis. In Sec. IIVI we move to study the 
strong coupling regime (U S> t a ) in order to determine a 
phase diagram which includes also different types of PS 
states. The results are summarized and the conclusions 
are drawn in Sec. W\ 



II. SINGLET-PAIRING TO CHARGE-DENSITY 
WAVE TRANSITION AT U < 

One of the open points raised in r ef. [l| is the exis- 
tence of regions at U < characterized by dominating 
charge-density wave (CDW) correlations instead of the 
singlet-superconducting (SS) ones that one has in the at- 
tractive symmetric HM. On a lattice the CDW and the 
SS correlation functions are defined respectively as 

C(r) = (rijUj+r) - (nj)(n j+r ) 



• \m > 



CDW 



Figure 1: Schematic representation of the phases discussed 
in the text. Large and small disks represent heavy and light 
fermions respectively. The ellipsis encircling the pairs in the 
SS phase denote onsite pairing. The HP state defined in Eq. 
p)l has not been drawn in real space because it is defined 
directly in momentum space by filling the two noninteracting 
bands up to /cf. 



P(r) 



(VjVj+r) 



where rij = + n 3 q while rjj — cj^Cff is the opera- 
tor that destroys an onsite pair with singlet spin wave- 
function. At large distances r, bosonization procedures 
predict Q 



C(r) 



K p cos(2kpr) 



A 



„K P +K a 



P(r) 



B 



„\IK„+K a 



where A and B are constants, and K p and K a are the 
Luttinger parameters for the charge and spin degrees of 
freedom respectively. Clearly C(r) dominates over P(r) 
when K p < 1, while K a has to be fixed to for gapped 
spin phases. A numerical estimate of K p from finite-size 
data can be obtained as in ref. [12] by considering the 
structure factor 



S(q)=J2j qr (^j+r) 



Here we have dropped the dependence on j because we 
implicitly assume that the correlation functions are trans- 
lationally invariant due to periodic boundary conditions 
(PBC). The value S(q = 0) corresponds to the average 
correlation and it diverges in the thermodynamic limit, 
since typically (rijiij+r) saturates at a finite value n 
for large distances. So one may consider the connected 
charge correlation C(r) in order to avoid this divergence. 
This choice affects the structure factor only at q = and 
bosonization predicts that the value of K p is directly re- 
lated to the limit q — > 

S(q) ~ K p — + ... 

7T 
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Following the procedure of ref. [l2| by selecting the 
smallest possible non-vanishing momentum compatible 
with PBC q± = 2tt/L one builds a sequence that approx- 
imates the linear slope 

Jf„(L) = S(gi)~. 

The dependence on L indicates that the sequence has 
to be extrapolated to L — » oo in order to obtain the 
limit q — > and the parameter K p . The DMRG results 
corresponding to various fillings n = 0.4, 0.6, 0.8 with 
z = 0.2, 0.5, 0.65, 0.8 and u = -0.001, -0.5, -3 are 
reported in table HI in the caption we have reported also 
the relevant features of our DMRG numerical calcula- 
tions. For a fixed value of u < we always find that 
the extrapolated K p decreases with increasing z. From 
this grid of points we can have an idea of the SS-CDW 
transition curve by locating the points at which K p = 1. 
This has been done interpolating the data with splines. 
An example of this procedure is given in Fig. [21 while 
the global results are plotted in Fig. [31 

In ref. [l| (Fig. 1 therein) the authors report a trian- 
gular region obtained by means of bosonization at u < 
where the dominant correlations are CDW or SS depend- 
ing on the filling. Here we observe that the shape of the 
transition line is indeed dependent on n: for u < —0.5 the 
separation line might have both a negative and a positive 
slope, depending on the value of n. For small (negative) 
u the curve has always a negative slope. Because of the 
uncertainty related to DMRG and finite-size effects, our 
estimate of the transition points is limited to u < — 10~ 3 ; 
moving closer to u = would produce values of K p es- 
sentially always equal to 1 within the numerical error, for 
all values of z. so we have not pushed our analysis and 
conclusions closer to u = 0. 

Finally we should mention that a direct inspection of 
the charge correlation functions in real and in Fourier 
space reveals that the only characteristic wavenumber is 
2kp = 7m, where typically the structure factor displays a 
peak for u < 0, but there is neither FFLO behavior - as 
expected since we have selected balanced species - nor a 
collapse (predicted at sufficiently large negative u [Hj]). 

III. WEAK-COUPLING LIMIT 

From a quantitative point of view, bosonization cannot 
be conclusive about the location of the HP to PS transi- 
tion at z = 1. This is due to the fact that, while a contin- 
uum limit approach is justified close to z = 0, in the FK 
limit even a small U can involve processes away from the 
Fermi surface, and so the requirement \U\ <C mm[tf,t±] 
strictly reduces the range of reliability of this approach [1| . 

Numerical data [22j indicate that in the highly asym- 
metric regime phase separation appears above a critical 
value of U which approaches zero in the low-density limit. 

When approaching the FK limit z — > 1, the kinetic 
energy of the lighter species becomes the dominant term 
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Figure 2: Example of cubic splines interpolation to locate the 
transition point and the associated error, using the DMRG 
data of table [I] The horizontal solid segment at K p = 1 
indicates the reported interval. 



of the Hamiltonian at weak coupling, U <C tf. For U > 0, 
the competition between the light particle kinetic energy 
and the repulsive on-site interaction can drive the system 
into an instability with respect to phase separation Q: 
there exists a critical value ups(z,n) above which light 
particles will occupy a large region of the system where 
n il = 0, creating an effective pressure which confines 
heavy particles into a small region with density close to 1; 
this effect is reminiscent of the segregated regime present 
in the FK model in the repulsive regime. 

Different numerical and analytical methods have been 
proposed in literature to identify this phase transition 
between the HP and the PS/TSS regimes [H H Hfj . In 
the following, we will (i) apply a second order perturba- 
tion theory approach, first introduced [27j to study the 
ground state of the weakly interacting symmetric HM, 
to compute the energy of a homogeneous phase ground 
state of the AHM and (ii) compare HP and PS or TSS 
energies in order to detect the line of quantum phase 
transition as a function of the original model parameters 
z, u, n. Finally, a comparison with previous results will 
be presented. 



A. Trial wave functions 

In the low coupling regime, we can consider as ex- 
tended ground state the exact one at U = (and ij. ^ 0), 
which can be obtained by filling both Fermi bands up to 
k Fa : 

i*>*p= n ii n 4 t i°) (2) 

\q\<k Fl |fe|<fe Ft 
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u \ z 




0.2 


0.5 


0.65 


0.8 




n = 0.4 










-1(T 3 




1.001 ±0.002* 


0.9978 ± 0.0002* 


0.994 ±0.001* 


0.986 ± 0.004* 


-0.5 




1.0700 ± 0.0004 


1.038 ±0.001 


0.992 ± 0.002 


0.906 ± 0.003 


-3 




1.4209 ± 0.0005 


1.2528 ± 0.0003 


1.0976 ± 0.0002 


0.86715 ± 0.00005 




n = 0.6 










-10~ 3 




1.0034 ± 0.0008* 


1.001 ±0.001* 


1.000 ± 0.002* 


0.998 ± 0.004* 


-0.5 




1.050 ±0.001* 


1.04 ±0.01* 


1.03 ±0.02* 


1.00 ±0.05* 


-3 




1.26016 ± 0.00007 


1.110 ±0.003 


0.969 ± 0.004 


0.751 ± 0.004 




n = 0.8 










-10~ 3 




1.008 ± 0.002* 


0.9976 ± 0.0001* 


0.983 ± 0.005 


0.951 ± 0.009 


-0.5 




1.045 ± 0.005* 


1.014 ±0.001 


0.980 ± 0.002 


0.913 ± 0.002 


-3 




1.1603 ± 0.0001 


1.0163 ± 0.0008 


0.880 ± 0.001 


0.6722 ± 0.0008 



Table I: Extrapolations for the parameter K p from DMRG simulations with PBC, L = 10,20,30,40,50, 1100-1300 optimized 
states and seven finite-system sweeps. This conservative choice guarantees an energy relative error of O(10 -6 ) up to L = 30 
and O(10~ J ) up to L = 50 (recall that the charge degrees of freedom are always gapless). Unless otherwise specified the 
extrapolations have been performed using quadratic fits in 1/L and the error bars are evaluated according to [2fJ (Chapter 15) 
using the sum of squared differences normalized to the fit degrees of freedom as a measure of the spread in the ordinates. %: 
Oscillating about reported value with spread. *: Linear fit in 1/L. 
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Figure 3: Numerical estimates of the transition lines from SS 
(on the left) to CDW (on the right) dominant correlations in 
the attractive regime for the values of filling reported in the 
legend. Not all the data used to obtain the figure have been 
reported in table U The error bars associated to the transi- 
tion points have been determined by means of splines passing 
through the upper and lower edges of each interval of possi- 
ble values of K p reported in the table (see the construction of 
Fig- O; when the uncertainty in z turned to be smaller than 
10~ 2 we have conventionally set it to this value to account for 
the approximation introduced by cubic splines interpolation. 
The lines joining the points are guides for the eye. 
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where |0) represents the zero-fermions vacuum and 
are the creation fermionic operators in Fourier space. 
While \iP)hp is n °t an exact eigenstate of the full Hamil 



tonian, perturbation theory above this ansatz have pro- 
vided excellent results for z = j2?J , where a comparison 
with the exact solution is possible, and, as shown later 
in the section, even in the highly asymmetric regime. 

The PS ground state instead can be obtained in the 
following way: first, we confine all heavy particles in 
a given part of the lattice of relative length v, with 
(L — Nf)/L > v > N^/L; then we consider two dif- 
ferent chains of length = l/L, = L(l — v) respec- 
tively and then fill the new Fermi bands till the momenta 
k' Fa = nN a /L a : 



!*(")> 



PS 



n 

\q'\< k Fi 



n 

\k'\<k' Ff 



fc't 



|0). 



In practice we have to consider that the effective light- 
particle and heavy-particle "chain lengths" are not L but 
L(l — v) and Li>, respectively. 

In addition, we can define a TSS as the one with a 
completely full region of heavy particles, i.e. v = N^/L. 
In this case, the variational wave function can be written 



a.s: 



I*) 



TSS 



(1- 



n 

ii)L<]<L 



i n 

w\<k' Ft 



5 f 



|0>. 



where now the down-fermion creation operators are taken 
in real space representation. Both TSS and PS state trial 
wave functions are eigenstate of the Hamiltonian up to a 
boundary term which we neglect in the following L — > oo 
limit. 
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B. Ground state energies 

The instability of a homogeneous ground state towards 
a phase separated one can be analyzed by computing the 
corresponding zero temperature energy: 

Ps(y\H\y)p S c HP {y\H\<S/) H p 

ops = r > c-HP 



L 



L 



and by comparing them to get the phase transition hy- 
persurface in parameter space described by: 

£ps(n^,n^,U,z) - £ H p(n^,n^,U,z) = 

A similar criterion can be applied to distinguish between 
TSS and PS state without segregation, as described later 
in the section. We remark that in the PS region, due 
to the fact that within the two subchains of length L^^ 
the up and down particles do not overlap, the interaction 
term provides at most a boundary contribution which can 
be neglected in the thermodynamic limit. We will come 
back to this point later. 

We can then compute the PS state energy density £p$ 
considering only the kinetic term contribution. For a 
general v, the result is 



£ps{v) = -t(l + z) 



2(1 



■ sin 



irn^ 
l-v 



(3) 



-t(l-z)— sinf^i') 

7T V V J 

As already stated in Q (in particular Sec. 3 therein), it 
is possible to fix the lowest energy state with respect to 
v searching for a minimum of ((3]) at fixed densities and 
hopping rates. The corresponding condition d v £ps = 
becomes 



(1 + z) sin 



TTTlf 



- — — COS 1 



optimal v = v*{z) in the range v S {n/2, 1 — n/2) sim- 
plifies to 



(5) 





-f 


-nn \ 


K 2(l-u) J 






7T71 \ 


.2(1-") ) 



with f(x) = sin a; — xcosi. The upper limit for v is 
formally obtained by compressing the up-fermions so that 
1 — v = n/2. When v = 1/2 the numerator of the right- 
hand side vanishes and z reaches the minimum possible 
value z = 0. When v = v m i„ = n/2 the right-hand side 
has the value 



z{n) 



-f{im/{2-n)) 
+ f(nn/(2-n))- 



This function decreases monotonically with the filling 
from z(0) = 1 to z(l) = 0. If z > z(n), the minimum 
energy is attained at the lower limit v* = n/2, corre- 
sponding to maximum compression of the down-fermions, 
independently of the value of z. Clearly, high densities 
favor a TSS state, since a large amount of light parti- 
cles produce a sufficient pressure in order to compress all 
heavy particles in a very small region. When the density 
is low enough, this condition cannot be satisfied, except 
at very large mass imbalances, and heavy particles still 
contribute to the kinetic energy of the system. Fig. [4] 
shows an example of such a construction for n = 1/2. 

The energy density of the HP state receives instead 
both kinetic and interaction contributions: £hp = £ t + 
£jj. The kinetic part is equivalent to the ground state 
energy density of the non-interacting case 

2 2 
£p = —til + z)— sin(7rn+) — <(1 — z)— sin(7rni ) (6) 

whereas the interacting part in the weak coupling limit 
can be expressed as a series in U applying a second order 
perturbation theory (27| : 



= (l-z) 



sin cos I i 



V v J 



V v J 



(4) 



If there exists a value of v* with 1 — n-j- > v* > nj,, 
which satisfies this condition, then the lowest energy 
state is \^(y*)) p$; otherwise, the minimum of £ps lies 
on the boundary v* = n^ and TSS is energetically more 
favorable. The boundary between these two regions is de- 
scribed by the condition d£ps /dv\ v = ni = 0; the solution 
of this equation provides a characteristic anisotropy co- 
efficient z for given n^^ that turns out to be independent 
of u. We expect that z represents a good estimate for the 
phase transition between TSS and v ^ nj, states even at 
intermediate couplings. In fact, we find good agreement 
with the values corresponding to the (almost) horizontal 
lines plotted in Fig. 3 of ref. [221 ) . 

As said before, here we consider only the balanced case 
n t = n i = f° r which the condition that yields the 



£u = jhp(^\Hu\^)hp 



~ ( HP {V\Hu Et 1 _ h Hu\*) 



HP 



0{U 3 ) 



where (...) com indicates a sum over connected diagrams 
and Hq denotes the unperturbed Hamiltonian. The first 
order contribution is obtained rewriting the number op- 
erators in momentum space, thus obtaining 



Hp(y\Hu\y)HP = LUn fni 



(7) 



while the second order contribution can be computed 
evaluating Goldstone diagrams: 



£ 



(2) 
HP 



U 2 #(k Ft ,k h a) 
i(27r) 3 (1 + z) 



(8) 
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Figure 4: Upper panel: Plot of the function z(ri). Lower 
panel: Optimal value of v as a function of z in the whole 
interval [0, 1] for filling n — 1/2. 



dq 
o sin <? 



dk' 



dkx 



sin(k + kpf) + a sin(fc' + kpi) 



Integrating the previous expression numerically, we can 
give a quantitative estimate of the phase boundary near 
the FK limit. Furthermore, by comparing £hp up to 
second order (plotted in Fig. [3]) with previous numerical 
results [22| . we have a good check that at half filling the 
variational ground state © is a correct description of the 
system even at finite U/t^ < 3. 



C. Phase boundaries 

The weak-coupling phase diagram is generally charac- 
terized by two types of phase transition: one between the 
HP state and the PS region, and the other one between 
different types of spatially separated states. 




0.5 1.0 1.5 2.0 2.5 3.0 



Figure 5: Energy density of the HP state at n = 1 for different 
asymmetry parameters: upper red line t^/tf = 0.1, lower blue 
line t±/t t = 0.3. 



Combining Eqs.©, ©, © and ©, the first men- 
tioned phase transition line is determined by the equa- 
tion: 



/7m\ 2 — n 

2sm (_j_ ( l +2) __ sm 



n 

T 



^(fcF t ,fcFt>a) 



(9) 



Table [H] shows how the correlated energy factor depends 
on Uf^ and z = (1 — o)/(l + a). In general, the larger 
the asymmetry z is (the smaller is a) the larger is i9. 
An illustrative plot at half-filling is presented in Fig. [5] 
Furthermore $ approaches in the low-density limit and 
grows with the filling. 



a\k F i = im/2 


tt/3 


tt/4 


tt/6 


0.01 


0.028 


0.022 


0.013 


0.05 


0.027 


0.021 


0.012 


0.1 


0.026 


0.02 


0.012 


0.15 


0.025 


0.019 


0.012 


0.2 


0.024 


0.018 


0.011 



Table II: Numerical values for ^(fc^t, &Ft; a ) I (2 71 ") 3 below half 
filling. Notice that the term (1 + z) in the denominator of Eq. 
([S]) has not been included in the definition of ■& so to have a 
more direct comparison with ref. [27) |. 

We will consider first the case when the density is 
medium-high: n-j- = = n/2 > 0.25, when a tran- 
sition from HP to TSS state should always take place, 
being v* > n/2 unfavored. By inspecting Eq. ([5]) it 
turns out that the ground state in the weak coupling 
limit is always homogeneous, even for large z (a < 0.2, 
i.e. z > 2/3). This fact agrees with previous numerical 
results [22| showing that the TSS phase is present only 
for U/t^ > 2.5. Even close to quarter filling, the phase 
transition in the FK limit is predicted at U/t^ = 1.15, 
which is beyond the weak coupling regime we are consid- 
ering in this section. 
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0.2 0.4 0.6 0.8 1.0' 



Figure 6: Second-order perturbative contribution at half fill- 
ing i9(7r/2, 7r/2, a)/(2n) for different values of a; solid line is 
a guide for the eye. 



Let us examine now the low density regime: = nr. = 
n/2 < 0.25, where the TSS might be the ground state 
only at very small a. Now we have to determine the 
transition point from HP to TSS state in the FK limit 
(a very small) as well as to explore the possibility of a 
transition to a PS state with v > for larger values 
of a. We will examine first the n = 1/3 case as an ex- 
ample. In the highly asymmetric regime (a < 0.01, i.e. 
z > 0.98), we again find a transition from a HP to a TSS 
state which happens for values of U(n = 1/3) /tf ~ 0.3. 
As explained in the previous section, using Eqs. (J3|) and 
(j4]) we can determine the phase transition line between 
TSS and a PS state with v ^ nj,. This transition occurs 
for ti(n = l/3)/tf ~ 0.025, a value which turns out to 
be independent of U since no contributions from the cor- 
relation energy are present. These transition points are 
in goo d accordance with the numerical results reported 
in [22j, and allow us to complete a general weak cou- 
pling phase diagram in the highly asymmetric regime for 
n = 1/3, which is shown in Fig. [71 At smaller densities, 
the critical value U(n)/t^ at which one finds the transi- 
tion from the HP to the TSS state becomes smaller. For 
example U(n = l/6)/t t - 0.06, while U(n = l/ 12 )/*t < 
10~ 3 . We notice that it is impossible, within our pertur- 
bative approach (see Eq. (JSJ), to find for this coefficient 
a value equal to zero: its value decreases as n becomes 
smaller but stays always finite, going to zero smoothly as 
n tends to zero. In addition, the TSS is always unstable 
with respect to a PS one, the transition now appear- 
ing at lower asymmetries (i[.(n = l/6)/tf ~ 2.5 x 10~ 3 , 
t ; (n = l/12)/t t ~ 2.7 x 10~ 4 ), whose values are still 
essentially insensitive to U. 



IV. STRONG-COUPLING LIMIT 

In this section we analyze the case of strong repul- 
sive onsite interaction between fermions, corresponding 
to t a <C U, a regime which is of particular interest for the 
experimental realization of the symmetric (z = 0) model 
in a cold atom system [T^, in which a Mott- insulator 
phase at half filling was found. One of the questions 
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Figure 7: Weak-coupling phase diagram with n — 1/3: white 
region represents the homogeneous state (HP), green regions 
are different phase separated states. The phase transition 
lines are computed using Eq. Q (dashed) and Eq. ©(solid). 
In this section we have adopted the parameters U/tf = u/(l + 
z) and a = t±/tf = (1 — z)/(l + z) instead of u and z to 
facilitate the comparison with previous results in the weak- 
coupling region. 



left open by bosonization is what happens to the HP-PS 
transition curve close to z = 0. From the phase diagrams 
in refs. [l[ and [l3[ it is not clear whether it approaches 
a finite value when z — > or, conversely, if there is a 
characteristic value of z ^ at which it diverges, as in- 
dicated also by some data on short sizes in ref. [23] (Sec. 
IV therein). To study this regime, we will construct an 
effective Hamiltonian that is able to describe the AHM 
when the parameter t can be considered as a small per- 
turbation with respect to U, for a generic filling, by using 
the method of the flow equations, developed by Wegner 
[IH and applied to the HM in [l||. The advantages of 
using such a technique are extensively described in [l9j . 
We only remark here that it yields a very general proce- 
dure which, in a recursive way, allows to find an effective 
Hamiltonian at any order of perturbation, for arbitrary 
fillings and geometries. 

We start by decomposing the fermionic Hilbert space of 
the model into the subspaces with exactly k fermionic 
pairs (double occupancies): W = ®^l^yHk- The projec- 
tors Pfe on these subspaces are defined via the generating 
function 



N 



J [1 — (1 — x)n^n iX \ . 



k=Q 



The kinetic energy term for the spin a fermions, T a 

J 



by 5k 



> ojg-ojcr, can be decomposed into three parts T a = 
Toa + T + \ a + T_i(j, which change the number of pairs 



'/?! 



0, +1, — 1. In other words: T„ 



8 



Sfclo Pk+mTcrPk- In these expressions, we have intro- 
duced the sum over < ij > which denotes nearest- 
neighbors sites i and j (with the couples ij and ji counted 
once each) since the procedure we are going to discuss is 
generalizable to any dimension. More explicitly: 

T a = E [ n i<r C h C i<r n 3<f + i 1 ~ n is)4a C l<y( l ~ n ]°) 

<ij> 

T +la = Y nisc\ a c 3<y {\ - n js ) 

<ij> 

T—\a = ^ ] (1 — n is)c[ a c j<y n jS 
<ij> 

It is not difficult to verify that [N d ,T ma ] = mT ma , 
where N d = Y2i n^n^, reflecting the transition from the 
Hilbert space T-Lk to T-Lk+m- To discuss higher-order in- 
teractions terms, it is useful to introduce products of 
hopping operators, T^... ak (m) = T rnirTl ...T mkrTk with 
the index vector m = (toi,7T12, ...,TOfc). It is found 
that the commutator of such an operator product with 
Nd involves the total weight of the product, M(m) = 

^k 

=1 



El 



and generally reads 



N d ,T£\.„ k { m ) 



M(m)Tif. ) .. (Tfe (m). We want now to find an effective 
Hamiltonian which does not mix the different Hilbert 
space sectors Tik, i.e. which conserves the total number 
of local pairs, thus suited to study physical properties at 
energy and temperature scales which are well below the 
Hubbard energy \U\. To do so, we consider a continu- 
ous unitary transformation which allows to remove in- 
teractions with non-vanishing overlap between different 
Hilbert space sectors. Thus, the transformed Hamilto- 
nian depends on a continuous flow parameter I: 

H(l) = -J2t a e (T (l) + UN d 

a 

where the generalized kinetic energy term Q a (l) contains 
all order interactions which are generated by the trans- 
formation: 



Here F^... ak denote suitable coupling functions that have 
to be determined by asking that the unitary transforma- 
tion cancels all terms that do not conserve the number of 
local pairs. The flow equations for these coupling func- 
tions follow from the equation of the flow for the Hamil- 
tonian Il8ll: 



k=l 



{m} 



(10) 



= E^E^)^. ) . 

fe=l {m} 



,(i;m)r« (m) 



Now, after imposing both the initial conditions 

(F^(0;m) = 1 and F^..<r k (0; m) = for k > 1) and 
the symmetries related to hermiticity and particle-hole 
transformation c\ a — > Ci a , which reverses the sign of the 
hopping term, one can recast the original flow equation 
(fTTj) in a recursive set of coupled nonlinear differential 
equations. From these equations it is easy to see that all 
the terms which connect different sectors of the Hilbert 
space vanish in the limit I — » oo. At the second order we 
find that the effective Hamiltonian reads 



-t 



t 2^ 

<ij> 



<ij> 



J [sfSf+SlSt+AiSfSj-wnj/A) 

<ij> 



-UN d (12) 



where S*' v ' z are the spin operators at site i, J = 2t^ti/U 
and A = (t^+t 2 l )/2t t t i > 1. At half filling, n = 1, we get, 
according to the Takahasi's theorem [20j | . that the terms 
corresponding to the odd orders of the expansion (in the 
our case the first two lines) vanish and we find the same 
Hamiltonian obtained in [2lfl representing a spin chain 
with an anisotropy term A > 1 along the z— axis (XXZ 
chain): the spin excitations are gapped and the spin- 
spin correlators decay exponentially with the distance. 
Also, in the limit tf = tj, (symmetric HM) the anisotropy 
term A becomes 1 and we find the well-known Heisenberg 
Hamiltonian (XXX chain), as it should be. 

It is well known (2l[ that for z > the system is in 
the Neel-like phase, with non-vanishing charge and spin 
gap and true long-range order. Here we are interested in 
examining the two limiting cases z = 1 (FK model) and 
3 = (Hubbard model) to study, more precisely, (?) the 
phase appearing in the Hubbard model when U — > oo, 
i.e. the so called spinless fermions phase (SF), where 
the orientation of the spins loses its relevance since the 
doubly occupied sites are strictly forbidden, and (ii) the 
state predicted in the FK model where the two fermionic 
species are demixed. 



dH{l) 
dl 



= [Va(l),H(l)} 



(11) 



which has been written here by using the (antihermitean) 
generator of the transformation 



T)a(l) 



U 



A. Spinless fermions 

The SF state of the Hubbard model in the limit 

U — > oo at filling n = N/L < 1 and equally populated 
species n-j- = n± = n/2 is rotationally invariant and in- 
variant under the up-down exchange. In this case, the 



expectation value of the hopping terms of the Hamilto- 
nian (TT2")) reads as follows 

sf (*i [-* J2(4« c i+ i « + h - c -) - tz X^sW^ + h - c -) 



2t 



+ tz ^2( c ]i c j+H + h.c.)]|*)sF = sin(vrn) 

3 

As for the J-terms, we can borrow directly its expression 
from (A3) of [H]: 



sf{^\Y. 3 njn j+1 \V)sF _ 2 sin 2 (7m) 



(13) 



Moreover, since I^sf is SU(2)-invariant, we can write 

sf(*\ gj S}S? +1 \*) SF = i sfW Ej s, ■ S j+1 \*)sf 
L 3 L 

(14) 

for a — x, y, z even if the J-part in the strong-coupling 
Hamiltonian is anisotropic, where now, from (A3) and 
(A4) of of Ull, we find 



SFWEiSi-Sj+lWsF (1 



In 2 

L V 4 
Collecting all the pieces it is easy to see that 
SF (y\H s - c \y) S F 2t 



sin 2 (nn) 



Ssi 



At 2 (l - z 2 ) 
U 



n 2 - 



L 



sin 2 (7m) 



= sin(7m)+ 

7T 



A 



- In 2 



As an example, in Fig. [8] we plot the local densities of 
fermions obtained numerically on a chain with L = 60 
and open boundary conditions (OBC), z = 0.1, n = 0.2 
and u — 100. The two species tend to occupy alternate 
regions but the fraction of double occupation is still sig- 
nificant. The comparison with the total density profile 
for spinless fermions at the same equivalent filling shows 
that the SF state is a good description of the ground state 
in this case. We have verified that the same happens if 
the filling is increased up to n — 0.9, the other param- 
eters being unchanged. On the contrary, if we still fix 
u = 100, n = 0.2 but increase the anisotropy to z = 0.3, 
appreciable differences in the density profiles start to ap- 
pear. 



B. Spatially separated states 

In the strong-coupling approach we can actually formu- 
late a slightly more general form of the PS state with re- 
spect to that of Sec. IIII Al which. however, leads formally 
to the same analytical expressions. Let us consider a se- 
quence of M contiguous intervals Ijf, Ij± {J = 1, . . . , M) 
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Figure 8: Spatial density profiles for an open chain with 
z = 0.1, filling n = 0.2 and u = 100. The continuous curve 
labelled SF is the local density for an equivalent chain of spin- 
less fermions (i.e. Eq. (|17p with £ = L and p = nL), showing 
a very good agreement with the spinful numerical data. 



and a state | v I'}ps in which the up and down spins are 
separated in the sense that there are no doubly occupied 
sites, Iif contains only up spins, contains only down 
spins, then again with up spins and so on. Let Lj a 
and Nj a , respectively, the number of sites and the num- 
ber of electrons in the interval We then make the 
further strong assumption that each up or down interval, 
irrespective of its length, is equally filled meaning that 
= Njf/Ljf and — Nj^/Lji do not depend on j. 
Now. if 



vL 



M 

i=i 



(l-v)L 



M 



are the total lengths associated with the motion of down 
and up spins we have Nf = Ej= ^Y?t = Ej=i "t-^it = 
n^L(l — v) and similarly N± = n^Lv and so, for equally 
populated species n-j- = n/[2(l — v)\ and = n/(2v). 

We first consider the thermodynamic limit in the 
case in which the interface points (which are 2M in 
number) do not contribute to the bulk energy density 
(limi_ >00 2M/L = 0) and, at the same time, each inter- 
val is extensively large ( lirri£,_>. 00 Lj a /L > 0), so that 
for every interval we can use the formula for the kinetic 
energy density of free fermions without worrying about 
finite-size and/or boundary effects: 



ps(*\T s _ c \V)ps 



L 



2U . ( n 
sm 7T— - 

7T V 2(1-1/) 



L 



2*4, . 




sm 




7T 
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— L (l — v) sin 7T— 

7T V ' V 2(1 - I/) 



2f 



x . / n 
— issm tt — 
TT V 2v 



Note that we do not necessarily require that the heavy 
(down) species is fully compressed, meaning v = n/2. 
The value for v will be determined variationally in order 
to give the smallest possible energy at a fixed z, exactly 
as done in Sec. MI B[ Let us now calculate the energy of 
such a state. 

As far as the J-term is concerned we first note that the 
transverse part is vanishing. In fact both SJ and Sj are 

composed by spin- flip terms like Cj a Cja but each interval 
contains spins of only one specie. Next, all the A-term 
can be rewritten as 



A ( Sj Sj+i 



■4 [("it 



n jl)( n 



■j'+it 




20 30 40 
Site index 



-{n ot +n ol ){{n J+n +n j+n )] = - — {n j - t n j+u +n j ±n j+n ). 

(15) 

When the expectation value on |\I/)ps is taken, the up 
and down parts factorize and there can be non- vanishing 
contributions only when j and j + 1 are at an interface 
between two intervals carrying opposite spins. If, as as- 
sumed above, the number of interface points does not 
grow as L we can neglect these contributions in the limit 
L — > 00. Therefore the energy density of the PS state 
reads 



Ps(^\H s _ c \y)ps 



2U. 

- — (l-,jsm^- 2(l 



2t 



I ■ ( n 
— vsin tt — 

TT V 2V 



As anticipated, this expression coincides with Eq. ([3]) in 
the balanced case = n± = n/2. In Fig. Owe present 
two examples of the spatial density profile for large u and 
intermediate/large z, from which the spatial separation 
of the two species can be clearly inferred. Note also that 
the light specie occupies regions with a non-vanishing 
fraction of empty sites and an oscillating density profile 
(rijf) is seen. Nonetheless the local density in the inter- 
vals occupied by the heavy fermions does not reach 1, so 
in these case we do not have a TSS as instead, for exam- 
ple, in Fig. 3 of ref. [HI valid for u = 20, z = 2/3, filling 
n = 0.8 on 40 sites (reproduced in our calculations but 
not shown here). 



C. PS': Extensive number of intervals M = aL 

In order to treat also the case in which the number of 
interfaces scales as a finite fraction a of the total number 
of sites, we will assume that all the intervals hosting up 
spins are equally long and equally filled: Lj-\ = t-\, Njf = 
Pf\/j = l,...,Mso that £-f = (1— v)ja, and = n/(2a); 



Figure 9: Spatial density profiles for open chains with z — 0.5 
and z = 0.9 (offset by +0.4 for clarity), filling n = 0.2 and 
u = 100. 



similarly for the intervals with down spins t\. = v /a and 
Pl = n/(2a) = p-f. Note that, for equally populated 
species, necessarily the finite number of electrons in each 
interval is the same for up or down spins, while the finite 
lengths are in general different. The energy density of 
this type of phase separated state will therefore have the 
form 

£ps> = Y( K t+ K -l)-^X^(("jt)("j+u) + ("ji)K-+it)) 

where Kf,i are the kinetic energies of p = p^ = Pi up 
or down fermions on intervals of length £^ with OBC, 
while the A-term comes from Eq. (|T5"j) and now can- 
not be neglected. The on-site terms (rij a ) also have to 
be evaluated in the same fashion and will be localized 
at the left or right end of the intervals (with equal val- 
ues). Let us denote them by r] a ; we have a contribution 

("rightt)(«loft|) + ("rightl) («leftf) = 2 VtVi for eac h of ^ e 

M pairs of up+down intervals so that 



£ 



PS' 



a(«t +Ki - 2JAr) t n i ). 



The calculation of the kinetic energy and of the surface 
density for an effective open chain of I sites with p free 
fermions is given in the appendix, leading to: 



£ps> = «*{(! + z)il- 



sin 


~7r(2p+l)" 


2(£ t + l) 


sin 




[2(^t+l). 



-(l-z){l- 



sin 






.2(^4 


-1)J 


sin 


7T 
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l + z 2 \ sin (^) J 1 sin (^rr) J 

(f t + l)(£; + l) 

The conditions p < £f and p < £± define the range of 
v £ [n/2, 1 — n/2], while the conditions p > 1,1^ > 1 
and £i > 1 imply a < min(n/2, i/, 1 — v ) = n/2. Once 
this expression is minimized by suitable values of a and 
i/ in this range we should, at least, check if the resulting 
energy density is smaller than the optimal energy density 
£ps determined above for the same values of z, n and, 
now, also u. 

Finally, we mention that we have also tried to enlarge 
the set of trial/ variational states by including the homo- 
geneous one (defined in Sec. IIIip . which is the correct 
ground state in the limit [7 = for all z. However we 
have verified that this additional state, for the fillings we 
have considered, could become relevant only when u < 1, 
outside the domain of validity of the strong-coupling ap- 
proach. Therefore, for the sake of compactness, we do 
not report these results here. 

D. Phase boundaries 




I i I i I i I i L_ 

0.2 0.4 0.6 0.8 



Z 

Figure 10: Transition lines in the (z, w)-plane at filling n — 0.4 
indicating how the ground state changes from SF to PS. In the 
separated regime above the transition lines, the state can be 
either PS (left) or TSS (right) (see Sec. |HTB| and the edge 
between the two is marked by a vertical line (analogously 
to Fig. 3 of ref. [12]). The dashed line corresponds to the 
bosonization prediction. 



From the condition £ps < £sf we get 

2(1 + 2) 



-(1 — v*) sin I 7r 



2(1-!/*) 



2(1 -z) 



- z) » . / n \ I . . . 
v sin 7r < sinlTm) 



4(1 - z 2 ) 



2^ r 



2 sin (im) 
n s 



2 + A / 1\ 1 

— - k2+ 4 -4 



that is u < ups(z) for d(n,z) > or u > ups(z) other- 
wise, having defined 



d(n, z) = sin(7rn) — (1 + z)(l — u*) sir 



2(1-1/*) 



— (1 — z)v* sin ( 7T 



and 



u PS (n,z) = 2tx{\-z 2 ) 



2 _ sin (7rn) 



(-1*2 + 1)- 



d(n, z) 



Now we can draw a phase diagram in the (z, u)-plane 
for a fixed value of the total filling n by indicating the 
regions where the PS or the SF state has the lower en- 
ergy. We show two examples (for n = 0.4 and 0.9) in 
Figs. [T0l and m~l respectively. We have analyzed in detail 
also the phase diagram for n = 0.6 (not shown), that 



turns out be qualitatively similar to the one for n = 0.4. 
Having in mind that bosonization could be considered 
quantitatively reliable only for small values of the inter- 
action, we have reported in the figures (dashed lines) the 
curves of Wentzel-Bardeen instability [13J where the ve- 
locity of one of the bosonization modes vanishes thereby 
indicating phase separation. In addition, we have consid- 
ered two typical cases of PS at n — 0.9 [L = 20,40,60), 
namely those at w = 100 for z = 0.1 and at u = 5 
for z = 0.9. In both cases, the charge structure fac- 
tor S(q) (as defined in Sec. |H]) displays a divergence 
for q — > typical of PS states [12| and a peak at 
q = 4kF (mod27r) = 27rn(mod27r) = 2tt(1 — n). 

In the strong-coupling approach the most interest- 
ing thing to understand seems to be the divergence of 
the transition line separating PS from SF behaviour 
at small z and large u. We have verified that when 
v = v*(z) is inserted into d(n,z) the denominator ap- 
pearing in ups(n, z) is always negative. In order to esti- 
mate ups(n,z) analytically we set v*{z) = 1/2 — zv\ — 
z 2 i*2+0(z i ) with v\ t 2 > (see Fig. [5]) and expand d(n, z) 
for z — > 

d(n, z) = -2i/i[/(7ro) - it 2 n 2 vi sin(7m)]z 2 + 0(z 3 ) 



so that at leading order 



Ups(n,z) 



1 



7T In 2 



[l + 0(z)] 



z 2 Ui(n)[f(Tm) — n 2 n 2 vi(n) sin(7rn)] 



(16) 



with f(x) = sin a; — a; cos a; as before. As far vi(n) 
as is concerned, by inserting v = 1/2 — zv\ + 0(z 2 ) 
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Figure 11: Same as Fig. [10] but with filling n = 0.9. Now 
also the PS' (see text) state is relevant for moderate values of 
u and the corresponding transition lines with SF and PS are 
marked. 

in Eq. ([5]) and solving at first order in z we get 
v\(n) = f(Tm)/[f'(irn)2Trn]. In summary, our analyti- 
cal approach predicts that there is no finite value of z 
below which PS disappears; by moving to a sufficiently 
large repulsive coupling it is always possible to induce a 
PS at arbitrarily small anisotropy. 



E. Inclusion of phase separated states with an 
infinite number of interfaces 

We have compared the energy density of the PS state 
at given values of n and z with the corresponding value 
for the PS' state discussed in Sec. IIV CI At small filling, 
say n = 0.1, the two variational solutions (with respect 
to v or v and a, respectively) coincide in the sense that 
the optimal value a* — > and the optimal value of v is 
the same. Moreover Bps' > £ps- When the filling is 
increased to, say, n = 0.4 the situation is similar with 
the exception of a small region {u < 0.8 for z = or 
u < 2.8 for z = 0.9) that can be considered to be outside 
the scope of the strong-coupling approach. At n = 0.6, 
the PS' solution can be ignored for u > 3.2 (a result 
checked at z = and z = 0.1) or u > 5.2 (as checked at 
z = 0.5 and z = 0.9). Thus, close to half-filling the PS' 
becomes relevant also at intermediate values of u and we 
have examined it in more detail. 

Let us fix n = 0.9; in the symmetric case z = the 
optimal value of a remains at its maximum a* = n/2 for 
u < 32 where it starts to decrease to reach a* = at 
u = 44.3; the PS' state has a lower energy density with 
respect to the PS one for u < 45.8. For positive z as long 
as z < 0.43 the PS' solution is never better than the ones 
considered before. When z increases further the PS' state 



is favored over the PS or even the SF one; the region at 
large z and moderate u where | v I'}ps' has a lower energy 
is characterized by the fact that a* = n/2 and v* = n/2 
so that p = £j_ = 1 meaning that all the down spins are 
isolated from each other. This configuration resembles 
the trimer crystal phase found in ref. [30] . with a mixture 
of hardcore bosons with attractive interaction and fillings 
1/3 and 2/3, which is equivalent to a repulsive case with 
balanced species and total filling 2/3 when a particle- hole 
transformation is performed. 

V. CONCLUSIONS 

Our study, which combines analytical calculations in 
the strong- and weak-coupling regimes and DMRG sim- 
ulations both for attractive and repulsive interaction, 
sheds some light on three qualitative and quantitative 
questions that are still open in the literature of the ID 
AHM: 

1. The shape of the transition line from SS to CDW 
dominant correlations for U < is filling-dependent 
and re-entrant in some cases (see Fig. [3]); 

2. Phase separation and phase segregation take place 
close to the Falicov-Kimball limit above an inter- 
action value U c which depends on the population 
in such a way that it approaches zero in the small 
density regime. Furthermore, transitions between 
phase separation and phase segregation at varying 
interaction take place at a nearly constant asym- 
metry; 

3. For small asymmetry, close to the Hubbard limit 
t-l ~ *T> the SF-PS transition takes place at larger 
and larger values of U ; Eq. ([To]) , obtained in the 
framework of a variational strong-coupling argu- 
ment, indicates that an arbitrarily small asymme- 
try is sufficient, at very large repulsions, to create 
a phase separated state which destroys the spinless 
fermion-like ground state of the Hubbard model. 
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Appendix: Free spinless fermions with open 
boundary conditions 

The eigenfunctions of the hopping operator 
— t X^=i 2 ( c j c j+i + h.c.) have the form 



fmU) = \ TTT sin ( fc ™i) i k ™ = -r—r , m =!,...,£ 
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and the dispersion relation is formally the same as in the 
case of PBC e(k m ) = — 2icos(fc m ) so 



COS 



rn — 1 



= t(i- 



7r(2p+l) 

2(i+i) 



2(^+1) 



where p < I is the number of particles. 

To compute the average density on the j-th site (rij) 
we pass to the creation/annihilation operators in fc-space 
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sin 



e + r 



The only non-vanishing possibility within the matrix el- 
ement for the vacuum |0) filled up to the momentum 
k p = up/ (I + 1) is m' = m. so we have the characteristic 
function of the Fermi sea nk m 
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The density at the edge is obtained by setting j = 1. 
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